# install.packages("xfun") #uncomment if need to install xfun
xfun::pkg_attach2(c("readxl", "dplyr","haven","stargazer","writexl","lmridge"))

load("Table01.RData")

# Note on the models relating to foriegn affairs ----
  # The survey uspew2004-10midpol, which has question relating to foreign affairs only,
    # was excluded from the models presented in the paper because of a mistake in the 
    # response rates provided by the pollsters in the original data.
  # objects relating to foreign affairs with an addendum .M indicates that this survey has been excluded from the data.
  # for robustness, we estimate the foreign affairs models a second time, including the survey with inaccurate response rates
    # these objects have the same name, without the addendum .M

# Note on dataset objects ----
  # objects starting with mydata are datasets used as input in the models.
  # the addendum .r indicates the data include a single unit response rate predictor.
  # the addendum .b indicates the data includes both contact and cooperation rates.
  # the following abbreviations indicate the topic that the data relate to:
    # econ = economy
    # civ = civil rights
    # energy = energy
    # imm = immigration
    # wel = welfare
    # foreign = foreign affairs

# Estimating the correct value of K ----

# the following code creates txt files with various calculations of K to be used in the ridge regressions
# The most appropriate calculation for our data is Muniz et al. 2009 (KM4)
# These estimations are very slowly and memory-demanding
  # we therefore have provided the txt files to allow users to view results without estimating 
  # users who wish to estimate the results again should uncomment the following lines (mark all, shift+command+c)
  # we recommend estimating on Macbook Pro with 32gb RAM. Computer with 16gb are likely to fail, especially the model on foreign data.

# mresponse refers to models with a single unit response rate predictor
# mboth refers to models where response rate is broken into two predictors: cooperation and contact.

# sink("keconomy20200806.txt")
# print("====================================================================================")
# mresponse.econK <- lmridge(abd_d ~ response + d_congress + year, data = mydata.econ.r, scaling = "sc", K = seq(0,10,0.001))
# print("abd_d ~ response + d_congress + year, data = mydata.econ.r")
# kest(mresponse.econK)
# print("====================================================================================")
# mboth.econK <- lmridge(abd_d ~ contact + cooperation + d_congress + year, data = mydata.econ.b, scaling = "sc", K = seq(0,10,0.001))
# print("abd_d ~ contact + cooperation + d_congress + year, data = mydata.econ.b")
# kest(mboth.econK)
# sink()
# 
# sink("kciv20200806.txt")
# print("====================================================================================")
# mresponse.civK <- lmridge(abd_d ~ response + d_congress + year, data = mydata.civ.r, scaling = "sc", K = seq(0,10,0.001))
# print("abd_d ~ response + d_congress + year, data = mydata.civ.r")
# kest(mresponse.civK)
# print("====================================================================================")
# mboth.civK <- lmridge(abd_d ~ contact + cooperation + d_congress + year, data = mydata.civ.b, scaling = "sc", K = seq(0,10,0.001))
# print("abd_d ~ contact + cooperation + d_congress + year, data = mydata.civ.b")
# kest(mboth.civK)
# sink()
# 
# sink("kenergy20200806.txt")
# print("====================================================================================")
# mresponse.energyK <- lmridge(abd_d ~ response + d_congress + year, data = mydata.energy.r, scaling = "sc", K = seq(0,10,0.001))
# print("abd_d ~ response + d_congress + year, data = mydata.energy.r")
# kest(mresponse.energyK)
# print("====================================================================================")
# mboth.energyK <- lmridge(abd_d ~ contact + cooperation + d_congress + year, data = mydata.energy.b, scaling = "sc", K = seq(0,10,0.001))
# print("abd_d ~ contact + cooperation + d_congress + year, data = mydata.energy.b")
# kest(mboth.energyK)
# sink()
# 
# sink("kimm20200806.txt")
# print("====================================================================================")
# mresponse.immK <- lmridge(abd_d ~ response + d_congress + year, data = mydata.imm.r, scaling = "sc", K = seq(0,10,0.001))
# print("abd_d ~ response + d_congress + year, data = mydata.imm.r")
# kest(mresponse.immK)
# print("====================================================================================")
# mboth.immK <- lmridge(abd_d ~ contact + cooperation + d_congress + year, data = mydata.imm.b, scaling = "sc", K = seq(0,10,0.001))
# print("abd_d ~ contact + cooperation + d_congress + year, data = mydata.imm.b")
# kest(mboth.immK)
# sink()
# 
# sink("kwel20200806.txt")
# print("====================================================================================")
# mresponse.welK <- lmridge(abd_d ~ response + d_congress + year, data = mydata.wel.r, scaling = "sc", K = seq(0,10,0.001))
# print("abd_d ~ response + d_congress + year, data = mydata.wel.r")
# kest(mresponse.welK)
# print("====================================================================================")
# mboth.welK <- lmridge(abd_d ~ contact + cooperation + d_congress + year, data = mydata.wel.b, scaling = "sc", K = seq(0,10,0.001))
# print("abd_d ~ contact + cooperation + d_congress + year, data = mydata.wel.b")
# kest(mboth.welK)
# sink()
# 
# sink("kforeign20200806.txt")
# print("====================================================================================")
# mresponse.foreignK <- lmridge(abd_d ~ response + d_congress + year, data = mydata.foreign.r, scaling = "sc", K = seq(0,10,0.001))
# print("abd_d ~ response + d_congress + year, data = mydata.foreign.r")
# kest(mresponse.foreignK)
# print("====================================================================================")
# mboth.foreignK <- lmridge(abd_d ~ contact + cooperation + d_congress + year, data = mydata.foreign.b, scaling = "sc", K = seq(0,10,0.001))
# print("abd_d ~ contact + cooperation + d_congress + year, data = mydata.foreign.b")
# kest(mboth.foreignK)
# sink()
# 
# sink("kforeignM20200806_response.txt")
# print("====================================================================================")
# mresponse.foreignK.m <- lmridge(abd_d ~ response + d_congress + year, data = mydata.foreign.r.m, scaling = "sc", K = seq(0,10,0.001))
# print("abd_d ~ response + d_congress + year, data = mydata.foreign.r.m")
# kest(mresponse.foreignK.m)
# sink()
# 
# sink("kforeignM20200806_both.txt")
# print("====================================================================================")
# mboth.foreignK.m <- lmridge(abd_d ~ contact + cooperation + d_congress + year, data = mydata.foreign.b.m, scaling = "sc", K = seq(0,10,0.001))
# print("abd_d ~ contact + cooperation + d_congress + year, data = mydata.foreign.b.m")
# kest(mboth.foreignK.m)
# sink()

# Estimating Ridge Regression ----

# the following code estimates the final models for each issue area using the correct value of K
# the models have already been estimated and appear as objects in the global environment
# users who wish to re-estiamte the models should uncomment the following lines

# unit response rate as single predictor:
# mresponse.econFINAL <- lmridge(abd_d ~ response + d_congress + year, data = mydata.econ.r, scaling = "sc", K = 3.31079)
# mresponse.civFINAL <- lmridge(abd_d ~ response + d_congress + year, data = mydata.civ.r, scaling = "sc", K = 2.32585)
# mresponse.energyFINAL <- lmridge(abd_d ~ response + d_congress + year, data = mydata.energy.r, scaling = "sc", K = 5.07074)
# mresponse.immigrationFINAL <- lmridge(abd_d ~ response + d_congress + year, data = mydata.imm.r, scaling = "sc", K = 4.71462)
# mresponse.welfareFINAL <- lmridge(abd_d ~ response + d_congress + year, data = mydata.wel.r, scaling = "sc", K = 3.04612)
# mresponse.foreignFINAL <- lmridge(abd_d ~ response + d_congress + year, data = mydata.foreign.r, scaling = "sc", K = 7.74354)
# mresponse.foreignFINAL.M <- lmridge(abd_d ~ response + d_congress + year, data = mydata.foreign.r.m, scaling = "sc", K = 8.00626)

# unit response rate broken into contact and cooperation:
# mboth.econFINAL <- lmridge(abd_d ~ contact + cooperation + d_congress + year, data = mydata.econ.b, scaling = "sc", K = 2.27576)
# mboth.civFINAL <- lmridge(abd_d ~ contact + cooperation + d_congress + year, data = mydata.civ.b, scaling = "sc", K = 2.82527)
# mboth.energyFINAL <- lmridge(abd_d ~ contact + cooperation + d_congress + year, data = mydata.energy.b, scaling = "sc", K = 3.92577)
# mboth.immigrationFINAL <- lmridge(abd_d ~ contact + cooperation + d_congress + year, data = mydata.imm.b, scaling = "sc", K = 4.41655)
# mboth.welfareFINAL <- lmridge(abd_d ~ contact + cooperation + d_congress + year, data = mydata.wel.b, scaling = "sc", K = 1.78701)
# mboth.foreignFINAL <- lmridge(abd_d ~ contact + cooperation + d_congress + year, data = mydata.foreign.b, scaling = "sc", K = 5.87731)
# mboth.foreignFINAL.M <- lmridge(abd_d ~ contact + cooperation + d_congress + year, data = mydata.foreign.b.m, scaling = "sc", K = 5.98371)

# Reviewing Results
summary(mresponse.econFINAL)
summary(mboth.econFINAL)

summary(mresponse.civFINAL)
summary(mboth.civFINAL)

summary(mresponse.energyFINAL)
summary(mboth.energyFINAL)

summary(mresponse.immigrationFINAL)
summary(mboth.immigrationFINAL)

summary(mresponse.welfareFINAL)
summary(mboth.welfareFINAL)

summary(mresponse.foreignFINAL.M)
summary(mboth.foreignFINAL.M)

# Robustness:
summary(mresponse.foreignFINAL)
summary(mboth.foreignFINAL)
